Social integration predicts survival in female white-faced capuchin monkeys

Abstract Across multiple species of social mammals, a growing number of studies have found that individual sociality is associated with survival. In long-lived species, like primates, lifespan is one of the main components of fitness. We used 18 years of data from the Lomas Barbudal Monkey Project to quantify social integration in 11 capuchin (Cebus capucinus) groups and tested whether female survivorship was associated with females’ tendencies to interact with three types of partners: (1) all group members, (2) adult females, and (3) adult males. We found strong evidence that females who engaged more with other females in affiliative interactions and foraged in close proximity experienced increased survivorship. We found some weak evidence that females might also benefit from engaging in more support in agonistic contexts with other females. These benefits were evident in models that account for the females’ rank and group size. Female interactions with all group members also increased survival, but the estimates of the effects were more uncertain. In interactions with adult males, only females who provided more grooming to males survived longer. The results presented here suggest that social integration may result in survival-related benefits. Females might enjoy these benefits through exchanging grooming for other currencies, such as coalitionary support or tolerance.


PATHWAYS
Below are some possible pathways through which grooming, foraging in close proximity of others, and coalitionary aggression interactions can result in survival related benefits. These pathways are based on pathways presented by Ostner and Schülke [1] and Thompson [2]. We chose the following survival-related benefits relevant to capuchin females: 1. Reduction of risk: less wounding during coalitionary aggression 2. Reduction of risk: less predation 3. Access to food resources 4. Social status: increase in rank Interactions in grooming, foraging, and coalitionary aggression can result in survival-related benefits as follows:

A female's interactions can directly result in survival related benefits
• Feeding tolerance results in access to food resources.
• Receiving coalitionary support results in less wounding and an increase in rank during rank-changing coalitions.
2. Female capuchins can exchange the same or different currencies with individuals in their social group • Grooming-giving exchanged for coalitionary support results in less wounding during coalitionary conflicts and/or an increase in rank during rank-changing coalitions [3,4].
• Grooming-giving exchanged for tolerance while feeding results in better access to resources [4,5].
• Grooming-giving in exchange for general tolerance results in less aggression and less wounding.
• Providing coalitionary support in exchange for tolerance while feeding results in more access to food resources.
• Providing tolerance while feeding in exchange for coalitionary support results in an increase in rank during rank-changing coalitions.
• Providing coalitionary support in exchange for coalitionary support results in less wounding, and/or an increase in rank during rank-changing coalitions. 3. A female's affiliative interactions might be a proxy for overall tolerance that the female enjoys from other group members. This can manifest as more tolerance in other fitnessrelated interactions.
• Feeding tolerance might be a proxy for the general tolerance of the female in other domains and might translate into less aggression and result in less wounding (interindividual tolerance in the context of food has been linked to social bonds [4,6] and increased levels of cooperation [7,8]. • Receiving grooming might be a proxy for the general tolerance of the female in other domains and translate into less aggression and less wounding. 4. A female's affiliative interactions might be a proxy for spatial centrality (the causal direction could go either way: interactions could cause centrality, or centrality might cause more interactions.) • Grooming (providing and receiving) is associated with spatial centrality and results in less predation.
• Feeding tolerance (providing and receiving) is associated with spatial centrality and results in less predation.
The pathways are not mutually exclusive and the list is not exhaustive. Finding an association between one interaction type (e.g. grooming-giving) and survival is consistent with multiple pathways.
The direction, i.e. being either a provider or receiver, matters in most pathways. For example, in exchanges of coalitionary support and grooming, it matters if the female is a recipient or a giver of coalitionary support. If she is a recipient of coalitionary support then she exchanges it by providing grooming. If she is a provider of coalitionary support, then she receives grooming in exchange. In the first case, the association between grooming-giving and survival would be consistent with the pathway to survival. In the second case, the association between coalitionary support-giving and survival would be consistent with the pathway to survival. Since direction matters in the pathways, whenever possible, we separated the individual roles of giver and receiver in order to understand how social interactions translate into increased survival.

CAUSES OF DEATH
For most females, we do not know the cause of death, because we rarely observe deaths as they occur. White-faced capuchin monkeys are female philopatric species. In our population we never observed a female disperse from the natal group, unless you count group fission events as a form of dispersal [9]. In a neighboring white-faced capuchin population at Santa Rosa, a handful of female dispersals were observed during a long-term study, but under very special demographic circumstances of alpha male takeover [10]. When a female vanishes, it is unlikely that she dispersed. We have demographic information surrounding these female disappearances, which further lowers the probability that female simply migrated completely out of our study area and was never seen again.
Out of 62 adult female deaths that occurred between 2002-2019, we have some information about the cause for 10 females (Table S1). Five females died during a drought and their deaths might be related to the drought conditions (possibly starvation). Two females were hit by a car, one was killed by a poacher, one disappeared when she was injured, and one vanished while pregnant. For all other females, we do not have any information about the cause of their deaths, because they usually disappear while the group is not being observed.
The absence of observations of infectious disease spreading through this population suggests that infectious disease is not a prominent cause of death for capuchins.
Old age is also not a common cause of death for females. The oldest living female in this population is estimated to be 39 years old, but the maximum life span in the wild is not known, because white-faced capuchins can live up to 55 years of age in captivity [11] and no field study has lasted this long. The distribution of deaths across age categories (Table S1) suggests that females die before they reach old age .
We have no evidence for differences in mortality patterns across social groups.

SOCIAL INTEGRATION ESTIMATES
A central challenge to estimating sociality is the heterogeneous sampling that characterizes ethological research [12]. Because some animals and their dyadic interactions are observed more than others, there is heterogeneous uncertainty about the underlying sociality of the individuals in the population. As an alternative to point estimates of sociality, our approach was designed to obtain measures of sociality that reflect that uncertainty, which subsequently carries forward into the accelerated failure time models.
To model social integration, we adapt the multilevel Social Relations Model, or SRM [13,14]. Originally formulated as an ANOVA model, the SRM partitions the variance of a directed dyadic outcome from actor i to partner j into (1) actor-level variation reflecting individuallevel heterogeneity of directing or sending ties to others, (2) partner-level variation reflecting individual-level heterogeneity of receiving ties from others, and (3) dyadic variation, which reflects the extent to which the i to j relationship is similar, or reciprocal, to the j to i relationship. The multilevel parameterization of the SRM advantageously allows this variation to be modeled with random effects, or varying intercepts, which accommodates imbalanced sampling relatively easily. Originally designed for continuous responses, the SRM can be adapted for binary, count, and ordered outcomes [15][16][17][18]. Implementations of these models have become common among researchers studying cooperative behaviour in humans [19][20][21].
In this case, the dyadic behaviour of the capuchin monkeys is a binomial outcome, where the numerator is the number of times that the behaviour was observed and the denominator is the maximum number of times it could have been observed given the sampling of the individuals and dyads. Using a "snapshot" approach that is common to behavioural ecology [22], we aggregated the observations by calendar year. Thus, an effect for each individual and dyad was estimated on an annual basis. 1 This timeframe was chosen primarily because the accelerated failure time models employed a similar annual timeframe. Each of the behaviours (grooming, coalitional support, and foraging proximity) was estimated separately, though we note that a multivariate SRM can be used to model multiple discrete behavioural outcomes [23].

A. Directed behaviours
Grooming and coalitional support are directed behaviours, where the behaviour originates from individual i and is directed toward individual j. The models assume that both individuals are members of group k in year t. These models can be notated: where α represents the intercept, a ikt is a node-level random intercept for the sending node i in group k in year t, b jkt is an analogous random intercept for incoming events to j in group k in year t, u |ijkt| is a symmetric random intercept for dyad ij in group k in year t, and d ijkt is a parameter that reflects the extent to which the directed i to j effect is distinguished from the symmetric effect. This latter parameter also serves to reflect possible overdispersion beyond the conditional expectation for the binomial distribution [24].
The random effects for sending and recipient nodes are assumed bivariate normally distributed with zero means and a homogeneous covariance matrix: As in other applications of the SRM [17], the correlation between the respective effects is known as the "generalized reciprocity correlation." In the present application, the modeling of this correlation helps to provide refined estimates of the node-level effects that are used as predictors in the survival analysis (i.e., the accelerated failure time models).
In practice, there is typically high correspondence between the node-level parameters and empirical calculations of out-degree and in-degree centrality [14]. That correspondence is evident in this study, and the primary benefit of using the Social Relations Model is to account for measurement uncertainty that stems from variation in the sampling effort of individuals. 2 The use of symmetric dyadic effects in this analysis parallels their use by Koster et al. [25]. That is, within each dyad, the relationship from ikt to jkt and the relationship from jkt to ikt share the same index, which is distinguished by vertical pipes |ijkt|. The effects are assumed to be normally distributed: u |ijkt| ∼ Normal(0, σ 2 u ) The effects for the directed dyadic effects are likewise assumed to be normally distributed: Note that binomial models lack the constant residual variance of Gaussian regression models. As an alternative, latent parameterizations of binomial models may assume that the corresponding variance is π 2 /3, or 3.29 [24].
The dyadic reciprocity correlation, ρ ijkt , can then be calculated by dividing the symmetric dyadic variance by the sum of the dyadic variances and the latent residual variance: This parameterization of ρ ijkt assumes that dyadic reciprocity must be positive, which we believe to be a reasonable assumption for the cooperative behaviours in this analysis. 3 In the context of the present analysis, the inclusion of the dyadic effects mitigates the possibility that node-level estimates are disproportionately reflective of, for instance, a high volume of interactions solely within dyadic contexts. Overall, the models show very high dyadic reciprocity correlations, which further mitigates concerns that the sociality measures from the node-level random effects might be largely a by-product of asymmetric interactions between individuals.

B. Symmetric behaviour
Unlike grooming and coalitional support, foraging in proximity to others is regarded as an undirected, or symmetric, behaviour. For symmetric behaviour, it is not appropriate to fit a Social Relations Model because the assignment of individuals to the respective i and j indices is arbitrary. The random effects for these individual should have a common variance, not the respective actor-level and partner-level variances that are estimated by the SRM [14,16,26,27]. 4 Individual-level variation may still be pronounced with symmetric behaviours, however, so we fit our model with individual-level random effects.
The binomial model examines the number of times that individuals i and j of group k were observed in proximity in year t, denoted y ijkt , as a proportion of the number of times that they could have been observed, n ijkt . The model is notated as: where α represents the intercept, v ikt is a node-level random intercept for the individual who has been assigned index i, and v jkt is a node-level random intercept for the individual who has been assigned index j. Note that the group effects are again assumed to be normally distributed with variance σ 2 m . A similar assumption applies to the individual-level effects, v ikt and v jkt . Importantly, note that they share the same variance, σ 2 v . Finally, an effect for the symmetric dyad, u |ijkt| , is included in the model to account for possible overdispersion in the variance of the binomial outcome [24]. As noted, we implemented these models in part to reflect the measurement uncertainty of heterogeneously sampled individuals. This uncertainty is evident in the individual-level estimates from the models, as seen in Figure S1. This figure, which depicts foraging proximity, represents the individuals from a single group in a single year of observation. Similar results are evident for other domains and other annualized observations of groups (not depicted here).  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Fig. S1. Comparison of the empirical data to model estimates for the domain of foraging proximity for a representative group of monkeys. Points correspond to individual monkeys who were present in the group during a calendar year. The top panel shows the empirical data, represented as the rate at which individuals were observed to be foraging in proximity to others in the group. The points are sized proportional to the number of observations of that individual in that year, per the legend. The bottom panel shows the model estimates for the same monkeys, as represented by the posterior means (i.e., the points) and standard deviations (i.e., the confidence bands, which represent 1 standard deviation). Note in particular that the model estimates show greater uncertainty (larger standard deviations) as the empirical sample size decreases for a given monkey. Owing to the partial pooling, the posterior means for these individuals have also been "shrunk" toward the population average.

DIRECTED ACYCLIC GRAPH (DAG)
Having many covariates in the model can result in "included variable bias" where predictors are not only causally influencing the outcome, but also influencing each other [29]. This can result in predictor induced statistical selection within the model and manifest itself through misleading statistical, but not causal, associations between the variables we are interested in. Our goal is to infer causal impact of social integration on survival. In addition, we considered (1) female's average dominance rank for that year, (2) her age, (3) average group size for that year, (4) the number of adult daughters co-resident with the female and (5) whether her mother was still alive.
An adult female's probability of dying in a given year is expected to increase with age [30,31]. Females who live in larger groups might experience greater survival rates due to reduced predation risks [32] and difference in interaction networks in comparison to females in smaller social groups due to a larger number of social partners available. Higher-ranking females possibly enjoy reduced mortality in comparison to lower-ranking females due to their central position in the group [33]. A female's rank and age are likely to influence how much the female participates in social interactions. Since white-faced capuchins are a female philopatric species, and females spend their entire lives with their female relatives, we considered including the number of adult daughters co-resident that year with the female and whether her mother was still alive [30]. The presence of a mother and/or adult daughters approximates a female's close kin network benefits: females who have more extensive kin networks might experience reduced mortality [31].
To check if the above assumed relationships justify the inclusion of all of the predictor variables in the model, we drew a directed acyclical graph (DAG) using the package daggity (v.3.0) and analyzed implied functional relationships. We found that in order to test for causal impact of social integration on survival, we need to include (1) female's average dominance rank for that year, (2) age, and (3) average group size for that year.

A. Age
Focal female ages were assigned based on demographic records of births and deaths collected from 1990 -2019. The ages of females who were born before 1990 were estimated in part by retroactively comparing photos taken then with photos of known-aged females collected later in the study. In addition, we inferred reproductive histories via genetic maternity data. We assumed that the age of first birth for each mother was six years, and that interbirth intervals were 2 years [9,34]. Out of 132 females in our sample, the birth year of 33 females was uncertain.

B. Average number of individuals in her group (group size)
This is the number of adult females, adult males, and immatures that resided in the female's group during the days when researchers spent at least six hours of observation with the group, averaged for the year. See Table S2 for group demography data.

C. Annual dominance index (dominance rank)
The annual dominance index represents the proportion of group members that the female dominated, on average, that year. For each observation day that the female (the focal) resided in a social group, we identified all of the other co-resident individuals (alters). To assess whether the female was dominant to an alter on a particular day, we first identified the dominance interaction immediately preceding and following the day of interest for each focal-alter dyad. Each focal-alter dyad could contribute 0, 1 or 2 data points to the daily dominance index. For each interaction, we identified that the focal individual was dominant if she was either the animal performing the supplanting, or being cowered at, or being avoided, or fled from. The daily dominance index, DDI i , of a focal individual, i, is a sum of dominance interactions where the focal was dominant to her alters, w i−a , divided by the total number of dominance interactions that the focal had with her alters, s i−a : Then, the average annual dominance index, ADI i , is an average of daily dominance indices: In some cases, either one or no dominance interactions were available for a focal-alter dyad. As a result, the individuals who did not have dominance interactions with the focal did not contribute to the calculation of the daily dominance index. Table S2.

AGE AND CHANGES IN INDIVIDUAL SOCIAL INTEGRATION
In this section, we provide plots that illustrate changes in social integration measures as a function of age. Across all dataset types (all group members, adult females only, adult males only) and social integration types (grooming giving, grooming receiving, support giving, support receiving, and foraging) the estimates tend to be lower for older females.

THE ACCELERATED FAILURE TIME MODEL
We specified the following model for the number of days before death, D i . The probability for the number of days before death comes from the cumulative probability distribution: For females who did not die during the observation period, the probability of surviving D i comes from the complementary cumulative probability distribution: We model the rate of dying, λ i , as follows: where α denotes the intercept or the base rate of number of days survived, α individual [id] denotes individual female random effects, and α group [groupid] denotes group random effects corresponding to the observation period. The model coefficients b sociality , b rank , b age , and b groupsize describe the impact of social integration, rank, age, and group size, respectively. We adopted a latent variable approach to model the social integration estimates given that the individuals' social integration estimates are not point estimates, but rather posterior distributions with means and standard deviations that reflect the measurement uncertainty.
We assumed a Normal (8, 0.5) prior for a base rate of survival, α, which places most of the prior mass between 0 to 20 years with the mean of 8 years and a long tail allowing more extreme values. For fixed effects, we assumed Normal (0, 1) priors. For individual-level random effects, α individual [id] and α group [groupid], we use a Normal (0, 1) prior.

ACCELERATED FAILURE MODEL RESULTS
Reported quantities are posterior means, standard deviation and 89% credible interval. Model diagnostic metrics provided are an estimate of the number of independent samples (n_eff) and an indicatorR (Rhat4) of the convergence of the Markov chains to the target distribution.
Each table indicates which social integration measure (grooming giving, grooming receiving, support giving, support receiving, foraging) was used as a predictor of survival, as well as which dataset (all group members, adult females only, adult males only) and how many adult females contributed to the survival analysis.

A. All group members dataset
This dataset includes adult female social integration measures estimated using their interactions with all group members.

B. Adult females only dataset
This dataset includes adult female social integration measures estimated using their interactions with adult females only.

C. Adult males only dataset
This dataset includes adult female social integration measures estimated using their interactions with adult males only.

ALTERNATIVE SPECIFICATIONS OF THE ACCELERATED FAILURE TIME MODELS
In this section, we present model results when we alter the parameterizations of the accelerated failure time models. We compared three types of models for the dataset containing adult female interactions of all group members: one with no random effects and two with individual-level random effect priors assumed as Normal (0, 0.2) or Normal (0, 1). All three model types resulted in similar estimates; therefore here we are reporting the results with individual-level random effect prior Normal (0, 1) in the main text, while the results of the other two types of models are reported here (Table S18  and Table S19). Table S18. Models without random effects. Each column represents the accelerated failure time model where social integration measure is either grooming giving, or grooming receiving, or support giving, or support receiving, or foraging. Estimates of fixed effects of each of the accelerated failure time models: posterior means and 89% HPDI.

Parameters
Grooming giving

ASSOCIATION BETWEEN PREVIOUS YEAR'S SOCIAL INTEGRATION ESTIMATE AND SURVIVAL
In this section, we present the model results when we model survival in year t as a function of the female's social integration measure in the previous year, t-1. The sample sizes of the datasets for these analyses are smaller than the datasets reported in the main text because in some cases behavioral data from the previous year were not available. Table S20 shows the results for the dataset containing female interactions with all group members. For this analysis, all groups contributed at least 6 calendar years of data (mean = 12.3, SD = 4.3, range: 6 -17). The subjects of this analysis were 125 females (899 female years), of whom 59 died during the observation period. Table S21 shows the results for the dataset containing female interactions with females only. For this analysis, all groups contributed at least 6 calendar years of data (mean = 11.5, SD = 4.4, range: 6 -17). The subjects of this analysis were 125 females (835 female years), of whom 55 died during the observation period. Table S22 shows the results for the dataset containing female interactions with males only. For this analysis, all groups contributed at least 6 calendar years of data (mean = 12.1, SD = 4.2, range: 6 -17). The subjects of this analysis were 125 females (882 female years), of whom 59 died during the observation period.
Overall, the inferences from this analysis are comparable to the inferences from the models reported in the main text. Table S20. Models with female's previous year's social integration measured through interactions with all group members. Each column represents the accelerated failure time model where social integration measure is either grooming giving, or grooming receiving, or support giving, or support receiving, or foraging. The reported quantities are posterior means (89% credible intervals in parentheses).

Parameters
Grooming giving